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Abstract 

How to distinguish and quantify deterministic and random influences 
on the statistics of turbulence data in meteorology cases is discussed 
from first principles. Liquid water path (LWP) changes in clouds, as 
retrieved from radio signals, upon different delay times, can be regarded 
as a stochastic Markov process. A detrended fluctuation analysis method 
indicates the existence of long range time correlations. The Fokker-Planck 
equation which models very precisely the LWP fluctuation empirical 
probability distributions, in particular, their non-Gaussian heavy tails is 
explicitly derived and written in terms of a drift and a diffusion coeffi- 
cient. Furthermore, Kramers- Moyal coefficients, as estimated from the 
empirical data, are found to be in good agreement with their first prin- 
ciple derivation. Finally, the equivalent Langevin equation is written for 
the LWP increments themselves. Thus rather than the existence of hi- 
erarchical structures, like an energy cascade process, strong correlations 
on different time scales, from small to large ones, are considered to be 
proven as intrinsic ingredients of such cloud evolutions. 



PACS numbers: 05.45.Tp, 05.45. Gg, 93.30.Fd, 89.69.+x; 02.50.Le, 05.40.-a, 
47.27.Ak, 87.23.Ge 



1 Introduction 

In order to establish a sound understanding for any scientific phenomenon, one 
has from numerical data to obtain laws which can be next derived theoretically. 
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Some difficulty arises in particular in nonlinear dynamical systems because of 
the problem to sort out noise from both chaos and deterministic components 
[Bhattacharya and Kanjilal, 2000; Provenzale et al., f992]. Moreover, algo- 
rithms for solving the inverse problem for nonlinear systems with underlying 
unknown dynamical processes are notoriously hardly reliable [Theiler, f992]. In 
fact, to extract dynamical equations for chaotic- like data is an enormous chal- 
lenge [Rowlands, 1992]. Practically one is often led to empirical relationships. 
This is the case of the meteorology field where there is a widely mixed set of 
various (sometimes) unknown influences, over different time and space scales. 
Works of interest for sorting out deterministic ingredients in chaotic systems, as 
in [Hilborn, 1994; Ott, 1993; Schreiber, 1999; Celluci, 1997; Davis ct al., 1996a], 
can be mentioned for general purpose though this is not an exhaustive list of 
references. 

It is known that two equivalent master equations govern the dynamics of a 
system, i.e. the Fokkcr-Planck equation and the Langevin equation, the former 
for the probability distribution function of time and space signal increments, 
the latter for the increment themselves [Reichl, 1980; Ernst et al., 1969; Risken, 
1984; Haggi and Thomas, 1982; Gardiner, 1983]. They are both condensates of 
the huge set of (6N) Hamilton equations which should in practice describe the 
whole dynamics of the system of N particles by giving the time evolution of both 
coordinates and momenta of each individual particle. This is a highly unrealistic 
scheme of work, and therefore conservation laws are used in order to derive the 
Navier-Stokes equations, from the averaging of basic quantities weighted by the 
above probability functions [Reichl, 1980; Huang, 1967]. Such a " micro-turned- 
macro-scopic description" can be by passed when describing the evolution of the 
quasi- equilibrium probability distribution function of such particles, i.e. in writ- 
ing a Boltzmann equation [Huang, 1967]. However this corresponds to a mean 
field description of the Fokker-Planck equation (FPE) [Reichl, 1980; Ernst et 
al., 1969; Risken, 1984; Haggi and Thomas, 1982; Gardiner, 1983], an approxi- 
mation which might be too rough for this sort of nonlinear system dynamics. 

Recently much advance has been made in describing nonlinear phenomena 
in meteorology, after the canonical Lorenz approximation [Lorenz, 1963; Schus- 
ter, 1984] of the Navier-Stokes equations for sorting out some deterministic and 
stochastic components of turbulent fluid motion in the atmosphere. In fact, 
turbulence [Frisch, 1995; Parisi and Frisch, 1985; Hunt, 1999] seems now to be 
part of more general class of problems similar to those found when there are 
long range fluctuations, as in financial data [Ghashghaie et al., 1996; Mantegna 
and Stanley, 1995], traffic [Chowdhury et al., 2000], dielectric breakdown [Van- 
dewalle et al., 1999], heart beats [Ivanov et al., 1999], etc... Such observations 
corroborate the idea that scaling [Stanley, 1971] and fractal geometry [Schuster, 
1984] principles could be useful in such studies on atmospheric turbulence [Mar- 
shak ct al., 1997]. Thus to derive the FPE [Risken, 1984; Haggi and Thomas, 
1982; Gardiner, 1983] for some meteorological effect is of primary interest. 

The aim of this paper is to derive such an FPE, in a first principle statistical 
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approach, from given raw data, i.e. for the Liquid Water Path (LWP) fluctua- 
tions in clouds, - in terms of a drift D^> and a diffusion coefficient. Much 
following the lines of thought of [Friedrich et al., 2000] it will be shown that 
the FPE (also known as the Kolmogorov equation) can be derived up to the 
first two moments of the (measured) conditional probability distribution func- 
tions p(Ax, At), where Ax is a signal increment and At a time delay. Thus the 
method is not anymore based on the conventional phenomenological compari- 
son between models and several stochastic aspects, but is a model independent 
(or first principle, from a statistical point of view) approach. This allows us 
to examine long and short time scales on the same footing, and suggests new 
features to be implemented in related weather forecasting. Indeed it is here 
below verified that the solution of the FPE yields the probability distributions 
with high accuracy, including the long and high tail Ax and At events. In some 
sense the agreement proves that searching for more improvement with including 
further event of similar cases is not necessary. 

Furthermore the so found analytical form of both drift and diffusion 
Z?( 2 ) coefficients has a simple physical interpretation. It is also shown that a 
truncation of the FPE expansion in terms of high order joint probabilities to the 
first two terms is quite sufficient. Moreover the identification of the underlying 
process leading to the heavy tailed probability density functions (pdf) for the 
LWP fluctuation correlation with At and the volatility clustering (as seen in 
Fig. 3 below) seems to be a new observation leading to an interestingly new set 
of puzzles for this sort of clouds [Cahalan et al., 1995], and more generally in 
atmospheric science. 

2 Data and theoretical analysis 

The changes in a time series for a signal x(t) are commonly measured by the 
quantity r — x(t + At)/x(t), or increments Ax — x(t + At) — x(t). Results of 
the analysis of a data set r(t), selecting among others LWP observations, are 
presented here. The 25 772 data points are taken from microwave radiome- 
ter measurements at the Southern Great Plains site of Atmospheric Radiation 
Measurement (ARM)[] of Department of Energy during the period January 9- 
14, 1998 as used by [Ivanova et al., 2000]. The microwave radiometer measures 
the brightness temperature at two frequency channels, one at 23.8 GHz and 
the other at 31.4 GHz. Then both brightness temperature data series are used 
to retrieve the vertical columnar amount of water vapor and vertical columnar 
amount of liquid water in the cloud, i.e. the so-called liquid water path. The 
raw data (i is equivalent to a time index) of the liquid water path in stratus 
clouds is shown in Fig. 1. 

The first thing is to consider whether the distribution function of increments 
has short or long tails, and whether there is some scaling law involved. Next 

1 http : 1 1 'www.arm.goi 
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the central issue is the understanding of the statistics of fluctuations which 
determine the evolution of the cloud. 

There are different estimators for the short and/or long range dependence of 
fluctuations correlations [Taqqu et al., 1995]. Through a detrended fluctuation 
analysis (DFA) method [Ausloos and Ivanova, 1999] we show first that the long 
range correlations are not brownian. The method has been used previously 
in the meteorological field [Ivanova et al., 2000; Ausloos and Ivanova, 1999; 
Koscielny-Bunde et al., 1998; Ivanova and Ausloos, 1999] and its concepts are 
not repeated here. We show on Fig. 2a the final result, for the 
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indicating a scaling law characterized by an exponent a — 0.36 ±0.01, markedly 
different from 0.5, thus indicating antipersistence of the signal from about 3 
to about 150 minutes. Next the power spectrum S(f) ~ .f _/3 with spectral 
exponent (3 — 1.57 ±0.02 is also shown (see Fig. 2b) [Davis et al., 1996b]. A 
Kolmogorov-Smirnov test on surrogate data has indicated the statistical validity 
of the value and error bars. 

Thereafter we focus on the LWP changes measured as increments Ax which 
are given in Fig. 3 in units of the standard deviation a of Ax at At = 640 s. 
In order to characterize the statistics of these LWP changes, LWP increments 
Axi, Ax2 for delay times At\, At 2 at the same time t are considered. The 
corresponding joint probability density functions are evaluated for various time 
delays Ati < At 2 < At 3 < ... directly from the given data set. One example 
of a contour plot of these functions is exhibited in Fig. 4. If two LWP changes, 
i.e. Axi and Ax 2 are statistically independent, the joint pdf should factorize 
into a product of two probability density functions: 

p(Axi, Atr, Ax 2 ,At 2 ) = p(Ax 1 ,At 1 ) P (Ax 2 ,At 2 ). (2) 

leading to a camelback or an isotropic single hill landscape. The tilted anisotropic 
form of the joint probability density (Fig. 4) clearly shows that such a factoriza- 
tion does not hold for small values of log(Aii/At 2 )|, whence both LWP changes 
are statistically dependent. The same is found to be true for other (Ati/ At j) 
ratios. This is in agreement with the hint taken from the above observations for 
the long range cross-correlation functions with the DFA. 

This implies a (fractal-like) hierarchy of time scales. If fluctuations in 
LWP went up over a certain At 2 , then it is more likely that, on a shorter 
Ati within the larger one, the LWP went down instead of up. To analyze 
these correlations in more detail, the question on what kind of statistical pro- 
cess underlies the LWP changes over a series of nested time delays Ati of 
decreasing duration should be raised. A complete characterization of the statis- 
tical properties of the data set in general requires the evaluation of joint pdf 's 
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p JV (Axi,Afi;...;Aa;jv,Atjv') depending on N variables (for arbitrarily large N). 
In the case of a Markov process (a process without memory) [Schuster, 1984], an 
important simplification arises: The N-point pdf p N is generated by a product of 
the conditional probabilities p( Axi+i, AU + i \Axi,Ati) = p(Axi+i,Ati + f,Axi,Ati)/ 
p(Axi,Ati) for i = 1,...,N-1. The conditional probability is given by the prob- 
ability of finding Ai i+1 values for fixed Axi. As a necessary condition, the 
Chapman-Kolmogorov equation 



p(Aa:i,At 1 |Ax2,Ata) = / d{A Xi )p{Axx, At^Ax^ Ati)p(A Xi , AU\Ax 2 , At 2 ) 

(3) 

should hold for any value of Ati, with Ati < Ati < At 2 . We checked the validity 
of the Chapman-Kolmogorov equation for different Ati triplets by comparing the 
directly evaluated conditional probability distributions p(Ax\, At\ \ Ax 2 , At 2 ) 
with the ones calculated (P ca i) according to Eq.(3), in Fig. 5. The solutions 
of Eq. (4) usually give rise to exponential laws, but initial conditions must be 
taken into account for adjusting the constants of the integration. Starting with a 
stretched exponential, the latter is conserved after successive integrations, with 
a width varying with At. The power law probability density function which is 
observed is thus a signature of fractal properties of the signal increments because 
the tails of the pdf are not truly exponential ones. In Fig. 5a, the contour lines 
of the two corresponding pdf's for all values of Axi are shown. There are mild 
deviations, probably resulting from a finite resolution of the statistics. Cuts for 
some exemplarily chosen values of Axi are shown in addition in Fig. 5 (b-d). 

As is well known, the Chapman-Kolmogorov equation yields an evolution 
equation for the change of the distribution functions p(Ax, Ai|Axi, Aii) and 
p(Ax,A£) across the scales At [Reichl, 1980; Ernst et al., 1969; Risken, 1984; 
Haggi and Thomas, 1982; Gardiner, 1983]. For the following it is convenient 
(and without loss of generality) to consider a normalized logarithmic time scale 
t= In (640/ At). The limiting case Ati ~ * corresponds to t — » oo. 

The Chapman-Kolmogorov equation formulated in differential form yields 
a master equation, which can take the form of a Fokker-Planck equation (for 
a detailed discussion, we refer the reader to [Reichl, 1980; Ernst et al.. 1969; 
Risken, 1984; Haggi and Thomas, 1982; Gardiner, 1983]f[ 



i p{Ax > T) = [ -sL° (1)(A ^ T) + i i)(2)(Aa:irMA ^ ) (4) 

in terms of a drift D^(Ax,t) and a diffusion coefficient D^ 2 \Ax,t). Note that 
the Fokker-Planck equation results from a truncation of the master equation 

2 The change in variable from At to r is non linear and might be debatable if the follow- 
ing equations (4) and (9) are integrated with respect to time. This is not done here. The 
integrations discussed in the text are made at fixed At or r. 
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expanded in terms of high order joint probabilities. According to Pawula's 
theorem [Risken, 1984], such a truncation is valid provided that the fourth order 
coefficient vanishes. We did check the coefficient for Ace G [—0.3,0.3] 
and obtained values that are two orders of magnitude smaller than and 
three decades smaller than D^-\ which we consider justifies the truncation. The 
functional dependence of the drift and diffusion coefficients can be estimated 
directly from the moments of the conditional probability distributions (cf. 
Fig. 5): 

M(k) = A7 / dAx '( Ax ' ~ &x) k p(Ax',T + At\Ax,t) (5) 
for different small At's ( Fig. 6), such that 

DW(Ax,t) = -J-limMW (6) 
k\ 

for At 0. 

The coefficient shows a linear dependence on Ax, while can be 
approximated by a polynomial of degree two in Ax. This behavior was found 
for all scales r and At. Therefore the drift term is well approximated 
by a linear function of Ax, whereas the diffusion term follows a function 
quadratic in Ax. For large values of Ax the statistics becomes poorer and the 
uncertainty increases. From a careful analysis of the data based on the functional 
dependences of and (Fig. 6 a-b), the following approximations hold 
true: 

D {1) = -0.0078Ax, (7) 

D {2) = 0.00259(Ax - 0.0305) 2 + 0.00009. (8) 

It may be worthwhile to remark that the observed quadratic dependence of 
the diffusion term is essential for the logarithmic scaling of the intermit- 
tency parameter in previous studies on turbulence. Finally, the FPE for the 
distribution function is known to be equivalent to a Langevin equation for the 
variable, i.e. Ax here, (within the Ito interpretation [Reichl, 1980; Ernst et al., 
1969; Risken, 1984; Haggi and Thomas, 1982; Gardiner, 1983]) 

^ Ax{t) = D {1) (Ax(t), t) + n(T)y[DM (Ax(t), t) , (9) 

where 7](t) is a fluctuating (5-correlated force with Gaussian statistics, i.e. < 
t](t) t]{t')> = 2S(t-t'). 
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3 Conclusions 



At least since the pioneering work of Lorenz [Lorenz, 1963] stochastic problems 
in turbulence are commonly treated as processes running in time t with long 
time correlations. Inspired by the idea of an existing energy cascade process 
[Cahalan, 1994] we present here a new approach, namely, we investigate among 
other related phenomena, how LWP changes are correlated on different time 
steps At. The pdf shape expresses an unexpected high probability (compared 
to a Gaussian pdf) of large LWP changes which is of utmost importance for 
forecasting. In recent works [Friedrich et al., 2000; Friedrich and Peinke, 1997] 
this finding leads to postulating the existence of hierarchical features, i.e. an 
energy cascade process from large to small time scales. The existence of finite 
non Gausssian tails for large events is thought to be due or to imply drastic 
evolutions, as for earthquakes, financial crashes or heart attacks. 

One unexpected result is the time dependence of the tails (Fig. 3) which seems 
quite smooth for the LWP case. This puzzle should initiate new research with a 
goal toward forecasting, in particular to find cases in which such a smoothness 
does not hold. This evolution shows how the pdf's deviate more and more from 
a Gaussian shape as At increases or r decreases. This definitely is a new quality 
in describing the hierarchical structure of such data sets, - not seen in the DFA 
nor spectral result. Now it becomes clear that one must not require stationary 
probability distributions for LWP differences for different r; on the contrary the 
coupling between different scales r via a Markov process is essential. Thus also a 
proper modeling of the time evolution of LWP differences for a fixed time delay 
must take into account the coupling of these quantities to LWP differences on 
different time delays. 

Furthermore, it is stressed that in contrast to the use of phenomenological 
fitting functions, the above method provides the evolution process of pdf's from 
small time delays to larger ones. Interestingly this is through an analogy with 
two physically meaningful coefficients, a drift term and a diffusion term 
Z)( 2 ) . The first one linearly behaves, thus looks like a "restoring force", the 
second behaving quadratically in Aa;, is obviously like an autocorrelation func- 
tion as for bona fide (chemical) diffusion. Of course further theoretical work is 
needed before understanding the numerical values in Eqs. (7)-(8). 

Finally, the present report presents a method on how to derive an underly- 
ing mathematical (statistical or model free) equation for a LWP cascade directly. 
The method yields an effective stochastic equation in the form of a FPE in the 
variable At. The excellent agreement between the experimental and the sta- 
tistical approach removes the need for much further proof based on examining 
many other cases for the same type of clouds. The FPE provides the complete 
knowledge as to how the statistics of LWP distribution change correlations on 
different delay times. Since this includes an autocorrelation analysis in time 
t for a scalar Ax, it is suggesting that the findings could be implemented in 
atmospheric weather low dimensional vector — models [Grabowski and Smo- 
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larkiewicz, 1999; Ragwitz and Kantz, 2000]. Consequently, even though for the 
first time in atmospheric science, a stochastic equation (of LWP evolution) is 
hereby introduced from first principles, several interesting statements can be 
presented and new questions opened. 
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Figure 1: Time dependence of a stratus cloud liquid water path (LWP) from 
Jan. 9 to 14, 1998, obtained at the ARM Southern Great Plains site with time 
resolution of 20 s. 
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Figure 2: (a) The detrended fluctuation analysis (DFA) function < F 2 (t) > x / 2 
(Eq. (1)) for the data in Fig. 1. Scaling range holds for a time interval between 
ca. 3 min and ca. 150 min. (b) Power spectrum S(f) for the data in Fig. 1. 
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Figure 3: Frequency p(Ax, At) of the LWP increments Ax for different time 
delays At; the units of Ax are taken as multiples of the standard deviation 
a = 0.0036 at At = 640 s. 
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Figure 4: Contour plot of the joint LWP increment pdf p(Ax\, At\; Ax2, A<2) 
for the simultaneous occurrence of LWP fluctuations Axi(Ati) and Ax2(At2); 
Ati — 80 s and At2 — 160 s. The contour lines correspond to logwp = 
-2,-2.5,-3,-3.5,-4. 
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Figure 5: (a) Empirical, and theoretical contour plots of the conditional LWP 
increment pdf p(Axi,Ati — Ax2,At2) for At2 = 180 s and At\ — 20 s. In order 
to verify the Chapman-Kolmogorov equation, the directly evaluated pdf (solid) 
is compared with the integrated pdf (dotted). Assuming a statistical error of 
the square root of the number of events of each bin, both pdfs arc found to 
be statistically identical; (b), (c), and (d) directly evaluated pdfs (dots) and 
results of the numerical integration of the Chapman-Kolmogorov equation (line) 
for cuts at Axi = -0.4, 0.0, 0.4. 
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Figure 6: Kramers-Moyal coefficients (a) D^> and (b) estimated from 

conditional pdf p(Ax 1} At 1 \Ax 2 , At 2 ); Aii = 20 s and At 2 = 40 s. The solid 
curves present a linear and a quadratic fit, respectively. 
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